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A CONSISTENT SPATIAL DIFFERENCING SCHEME FOR THE TRANSONIC 


FULL-POTENTIAL EQUATION IN THREE DIMENSIONS 
Scott D. Thomas* and Terry L. Holst 
Ames Research Center 


SUMMARY 


A full-potential, steady, transonic wing, flow solver has been modified so that 
free-stream density and residual are captured m regions of constant velocity. The 
numerically precise free-stream consistency is obtained by slightly altering the 
differencing scheme without affecting the implicit solution algorithm. The changes 
chiefly affect the fifteen metrics per grid point, which are computed once and 
stored. With this new method, the outer boundary condition is captured accurately, 
and the smoothness of the solution is especially improved near regions of grid 
discontinuity. 


INTRODUCTION 


Free-stream consistency is defined as the ability to numerically capture the 
free-stream density and vanishing residual from the transformed governing equations 
on an arbitrary computational mesh. The calculation of these quantities must be 
algebraically precise; that is, accurate to within machine roundoff, instead of to 
within some order of magnitude associated with the mesh spacing. 

Following the work of Flores et al. (ref. 1), modifications have been made to 
the three-dimensional (3-D) full-potential flow solver described by Thomas and Holst 
(ref. 2) in order to implement a free-stream-consistent spatial differencing 
scheme. This program uses a coordinate transformation from physical to computa- 
tional space in order to represent accurately the aerodynamic surface. It models 
the inviscid, compressible, irrotational flow of air past a lifting wing mounted on 
a wall. The vectorized computer code TWING, Cycle 3, has been upgraded to 
Cycle 4. It is designed to run on a Cray supercomputer. 

The spirit of the present approach follows the description in reference 1 in 
that it only requires careful numerical definition and use of metric terms required 
by the coordinate transformation. This is distinct from other approaches involving 
the subtraction of "artificial flux" terms as mentioned, for example, by Caughey and 
Jameson (ref. 3). 


*Informatics General Corporation, Palo Alto, CA 94303. 
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Beneficial features of TWING, Cycle 4, include: 

1 . Elimination of grid-induced perturbations in the flow field near the outer 
(or free-stream) boundary 

2. Reduction of perturbations in the flow field that otherwise would be intro- 
duced by grid discontinuities 

3. Improved ability to treat wings with sharp leading edges 

On the negative side, instead of seven metric quantities per point, now fifteen 
are required. However, a new data compaction and expansion technique keeps the 
storage requirement close to that for the previous cycle. 

The operation count per iteration per grid point remains about the same as 
before, reduced by the elimination of some averaging operations, but this is offset 
by the addition of the data expansion operations. The number of iterations required 
for convergence, which is usually taken as a two-order drop in the maximum modulus 
of the residual, increases for some cases, yet decreases for others. The initial 
residual for Cycle 4 is usually less than two orders of magnitude smaller than 
before. 


GOVERNING EQUATIONS 


This description is taken from reference 2 and is included here for future 
reference. The 3-D full-potential equation written in strong conservation-law form 
in Cartesian coordinates is given by 


(p<t> ) + (p 4> ) + ( p4> ) = 0 
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where the density p is scaled by the stagnation density p g , the velocity compo- 
nents 4> x , 4>y, 4> z are scaled by the critical speed of sound a*, and y is the 
ratio of specific heats. The value 1.4 is assigned to y in all cases under 
discussion . 


The orthogonal physical coordinates x, y, z are oriented in the flow, span, 
and vertical directions, respectively. A general coordinate transformation is 
applied to the above equations, and the implicit AF2 relaxation takes place in the 
new coordinate system. Equation (1), after the transformation, is 


(pU/J) + (pV/J) + (pW/J) = 0 
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where £, n, and z represent wraparound, spanwise, and radial-like directions, 
respectively, and the contravariant velocity components U, V, and W are defined by 
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The quantities H and J are defined by 


. a(g.n.c) 
a(x,y,z) 



(3) 


(4) 


and J = det(H). Dirichlet boundary conditions for the potential are imposed at the 
outer free-stream boundaries, and Neumann boundary conditions that yield flow tan- 
gency are imposed at the symmetry wall and wing surface boundaries. 

The discretization of the coordinate space replaces the real variables £, n, 
and z with integer variables i, j, and k, respectively. 


CONSISTENT FREE-STREAM DENSITY 


Assume that the free-stream potential function is given by 

4> = ux + vy + wz (5) 

CD CO 03 

T 

so that the gradient, v<)> = (u ,v ,w ) , gives the free-stream velocity vector. From 

OO Cu 00 

equation (1b) it is seen that the density is a function of the magnitude of the 
velocity vector, so that 


= p(u + v £ 


wf) 


( 6 ) 


The method used to capture the free-stream density is simply to difference <j> in 
i, j, and k the same way that x, y, and z are differenced in order to find the 
metrics. The metrics are defined to be the elements of the matrix H^H (see 
eqs. (3) and (4) above). An algebraic proof that this yields free-stream density is 
given below. This will also serve as an introduction to the next section on consis- 
tent residual computation. 
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From equation (2b), one sees that it will suffice to show that 

U4> + V<t> + W* = u 2 + v 2 + w 2 (7) 

5 n t <= => m ' 

when finite differences replace derivatives. Equation (3) can be used to establish 
that the inner product 

U<t> + V* + W<b = (U,V,W)(4> ,♦ ,♦ ) T 

s n v* s n v* 

= ( <J> r 9 1 4> ) H^*H( 4> , <J> ,4> 

s M s M 

= a t a (8) 

T 

where the column vector A = H(<t>_,4> , <t> ) . The method used to compute H is first 
to form a matrix of finite differences of x, y, and z, 


3(x,y,z) 
3(5, n,s) 


x 
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(9) 


and then invert to find H = K _1 . For instance, the quantity x^ at (i+1/2,j,k) is 
found to be x i+1fJ>k - x lJ|k . 

By the chain rule, the following holds analytically: 



( 10 ) 


When differenced numerically, this becomes an approximation. The symbol K* shall 
denote the matrix on the right-hand of equation (10) when the difference formulas 
for <)> are identical to those used for x, y, and z. Thus K* is the same as K 
of equation (9). Finally, it is seen that equation (3) and the definition of <t> 
leads to the identity 


A = K'Vu ,4) ,4) ) T = (u , v , w ) T (11) 

X 7 yZ CO 7 CD* CD ' 1 

Equation (7) follows from (8) and (11), which completes the proof. 

This technique has been implemented in TWING Cycle 4 to find the densities at 

the (i+1/2,j,k) grid locations. Six density metrics are required per grid point 

T 

since the metric array, H H, is symmetric. Whenever the density is needed at other 
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locations, like (i,j+1/2,k), the value is obtained from an average of neighboring 
points. In a free-stream region of flow all of the densities will be the same, so 
that the generality of the formulation is retained. 


CONSISTENT FREE-STREAM RESIDUAL 


This section examines how the metrics must be defined in order to have the 
numerically obtained residual vanish identically, independent of the shape of the 
grid. Since the Jacobian of the transformation appears in the denominator, we make 
the assumption that it remains positive. This corresponds to the condition that 
grid lines of the same family do not cross. The density is constant, so it may be 
removed from equation (2a), leaving: 


(U/J)_ + (V/J) + (W/J) = 0 (12) 

5 n C 

Treating this analytically, the left-hand side of equation (12) may be rewritten as 
an inner product: 


[det(K 1 )(K 2 


Once again, the letter K denotes a transformation matrix as defined in equa- 
tion ( 10) . The superscript and subscripts on K will be used below to indicate 
slightly different differencing methods used to approximate the analytically 
defined K. As in the previous section, K* denotes the matrix that arises from the 
method used to difference <t>. 

In equation (13) we actually have control over the part of the right-hand side 
that appears in square brackets. It is evident that by farming Kg with the same 
difference formulas as those used in forming <t>, then K” K* = I may be dropped 
from the equation. The free-stream condition on 4> then allows us to rewrite this 
expression as: 



Vk^Ik# 


(13) 
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It only remains to show how to form the remaining K matrices so that this expres- 
sion will vanish identically when finite differences replace the derivatives. 
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The values assigned to u D , v m , and w m are arbitrary, and equation (14) is 
linear, so it will suffice to treat the case for which u m = 1 and v = w =0. 

7 00 co co 

This amounts to an assumption about the direction of the free-stream velocity. The 
general case will follow by permutation. 

Make the assumption that = K 2 = K (drop the subscripts), since analytically 


(K-V . 


1 


det(K) 


y 2 - y z 

J n c n 

y c z 5 ' y 5 z c 

Vn - y n z C 


(15) 


so the determinant will divide out. Six elements of the matrix in equation (15) 
have been omitted because the assumption on the free-stream direction renders them 
immaterial. The problem has been reduced to showing how the following expression 
can be made to vanish: 
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(16) 


This is known as the "third consistency condition" in reference 1. It should be 
clear that, if this were an analytic expression, then it would indeed vanish by the 
property of equality of mixed partial derivatives. 


The first summand of equation (16), (y z - z y ) , is approximated by 

5 S 


( V< ' Vc’i.i/a.j.k - ( Vc - Vc’i-i/a.j.ic 


(17) 


A natural way to expand the first derivative in this expression, y^, is 
^i+1/2,j+1/2,k - yi+1/2,j-1/2,k)- Now replace y i+1/2f j+ i/ 2 ,k with an average 

of y at cell centers, namely (y i+1 /2, j+1/2,k+1/2 + ^1+1/2, j + 1/2,k-1/2 )/2 ' Such 
averages may be applied to all the other terms in equation (17). This procedure may 
be carried to the point where every term in equation (16) is replaced by finite 
difference expressions involving only the values of x, y, and z at the eight 
neighboring cell centers surrounding the grid point at (i,j,k). 


Complete cancellation can be obtained by forming the primitive quantities for 
x, y, and z at the cell centers using averages of the eight nearest neighbors. For 
instance, 


x i+1/2,J+1/2,k+1/2 


< 1 /®> < K i , j , k + i£ 1 +1,j,k + x i,j+l,k + x i+1 , j+1 ,k 
* x i,j,k+1 * x i-.1,J,k+1 * x i, J+1 ,k+1 + x U1,j*1,k+1 ) 


( 18 ) 
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As a point of interest, Pulliam and Steger (ref. 4) suggested a similar averaging 
for the Euler equations, but they studied this "third consistency condition" without 
the other features mentioned above. 

Verification of the cancellation property is left as an exercise for the 
reader. Numerical verification has been conducted with a small Fortran program 
called RFTEST. It was set up to run through a sequence of random perturbations of a 
39-point difference molecule surrounding the single point of interest. The resid- 
ual, the left-hand side of equation (2a), was found to be consistently close to zero 
when three different machine precisions were used. Four, eight, and sixteen-byte 
real-number representations were tested on a DEC VAX computer. 

It is customary to discretize equation (2a) as a sum of three flux terms (shown 
one per line) 


<‘ >U/J >i.1/2,J, I c - <<’ U/J) i-1/2,J,k - <*™>i,j + 1/2,k - (oV/J) i,J-1/2,k 


<» H/J) i,j,M/2 - '» W/J >l,j, t -V2 = 0 <'9> 


The flux terms at adjacent points are shared, so at each point of the mesh there are 
three "flux metrics" required for each of the three coordinate directions. Although 
three of these metrics are computed at the same half-points as the density metrics, 
they must be kept separate for free-stream consistency. Thus, six plus nine (or 
fifteen) metrics per point are computed in order to impose free-stream consistency 
for the full-potential equation in three dimensions. 


It is an interesting exercise to show that if the 3-D grid is composed of a 
stack of parallel 2-D grids, then some of the density and flux metrics may be shared 
while retaining free-stream consistency. At the same time, the compactness of the 
differencing molecule for density and residual may remain as described above. In 
this case, thirteen metric quantities would be needed per grid point. An early 
version of the RFTEST program was used to show that this property does not hold for 
a general computational mesh. This method was also tested in TWING but was later 
replaced by the fifteen-metric method described above. It should be noted that a 
remnant of the earlier method remains implemented in Cycle 4 m the wing extension 
region beyond the tip for the innermost k plane. 


In reference 1 it is shown that a consistent differencing scheme can be formu- 
lated with only thirteen metrics by forming the <j> differences for density in a 
manner consistent with the way the metrics are formed for the residual. This would 
lead to a larger density differencing molecule and, thus, to more numerical smearing 
of the solution in the vicinity of a shock; therefore, it is not recommended. 


PERFORMANCE NEAR A GRID DISCONTINUITY 


There are several cases in which it is convenient to use a computational grid 
that does not have optimal smoothness. A simple example of this is in the treatment 
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of a wing with break discontinuities in which the sweep angle of the leading and 
trailing edges may change abruptly. The consistency method guarantees good results 
using a simple grid which may propagate the break discontinuity to the outer bound- 
ary. In another example, one may wish to use a block approach m which a transonic 
wing flow solver just treats wing-like elements of a complicated geometry. In this 
case, it may be necessary to make the outer boundary of the wing block conform to 
the surfaces presented by the other blocks. 

Figure 1 shows an example of a grid about a wing that is mounted on a fuselage- 
like bump on a wall as an example of a case for which one might want to use a grid 
that is not made of a stack of two-dimensional grids. In this case, we may choose 
to propagate the bump throughout the grid from the root wall to the opposite free- 
stream boundary. The approach presented above applies in a limited sense to regions 
of the flow field where the velocity changes slowly, such as fore and aft of the 
wing. The region near the wing enjoys denser packing of grid points owing to the 
topology chosen for the flow solver, so consistent differencing is not as major an 
issue. 



Figure 1.- Grid about a wing on a fuselage. 


In the following discussion, computed results of Cycle 3 and Cycle 4 of the 
transonic wing analysis code TWING are compared. Cycle 4 of TWING uses the consis- 
tent spatial differencing described in the previous section, but the earlier Cycle 3 
version does not. 
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Comparison of Cycle 3 and Cycle 4 


The grid generator supplied with TWING is designed to build a three-dimensional 
grid by stacking two-dimensional grids, located at j = constant values, along the 
spanwise direction. Near the tip of the wing there is an abrupt transition from 
airfoil cross sections of positive thickness to zero thickness. The grid propagates 
this abrupt change from the inner to the outer boundaries, and its effect can be 
clearly seen by observing the velocity vectors in a fixed k plane. Figure 2 was 
generated by Cycle 3 with the default M6 wing at Mach 0.84 and zero angle of attack 
on a grid of 40,050 points (89 k 25 x 18). It shows a plot of the projections of 
the velocity vectors in the computational plane for k = 6 (about one root chord 
downstream) into a plane that is located aft of the wing that is aligned perpendicu- 
lar to the free-stream direction. Since the angle of attack is zero, we expect that 
the projections will have almost no magnitude. Observe that this does not happen at 
points lying in the three j planes just beyond the tip. In addition, the direc- 
tions of the projections change abruptly as we move from j planes that contain 
wing sections to j planes beyond the wing tip, which do not. In figure 3, we see 
the same case as above except that the velocity vectors were generated with 
Cycle 4. We see that each projection does have almost zero magnitude and also that 
the directions of the projections transition smoothly throughout the entire region 
that is shown. 



Figure 2.- Velocity vector projection for M m = 0.84, a = 0°, TWING Cycle 3. 


9 




Figure 3.- Velocity vector projection for = 0.84, a = 0°, TWING Cycle 4. 

Again we compare the two cycles of TWING in the next two figures, but in this 
case the angle of attack assumes its default value of 3.06°. Figure 4 shows the 
projections of the velocity vectors for the same k plane as before, using 
Cycle 3. Since the angle of attack is nonzero, the plane of projection has been 
inclined 3.06° so the wake region appears to be offset from the wing location. A 
pattern that resembles a vortex is generated at the edge of the wake region near the 
tip, despite the inherent potential assumption of flow irrotationality . The appear- 
ance of this flow pattern can be explained by the fact that a circulation model is 
added to the governing equations m order to simulate a lifting wing. This is a 
common approach among wing codes (flow solvers) that are based on a potential formu- 
lation. Upon comparison with figure 5, which was generated by Cycle 4, we see that 
the irregularities in the flow field near the tip have been eliminated. 


CONSIDERATIONS FOR A SHARP LEADING EDGE 


Some applications require the analysis of a wing with a sharp leading edge. 

This is considered to be another type of grid discontinuity, especially with the "0" 
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Figure 4.- Velocity vector projection for M ro = 0.84, a = 3.06°, 

TWING Cycle 3. 

mesh topology used in TWING. It was found that Cycle 4 (with consistent spatial 
differencing) is capable of treating a wing with a sharp leading edge, as opposed to 
Cycle 3 which is not. To demonstrate this point, a simple wing shape was devised 
having rectangular planform and an airfoil cross section for which the upper part is 
taken from a 12#-thick circular-arc airfoil and the lower part is taken to be a flat 
plate. This results in an airfoil that is 6% thick. 

Figure 6 shows a plot of the coefficient of pressure for two wing stations 
(root and mid-span) which illustrates that a solution generated by Cycle 3 exhibits 
erratic behavior near the leading edge discontinuity. This excessive gradient can 
cause stability problems under some conditions. The Mach number for this case is 
subsonic, at 0.2, and the angle of attack is zero. The leading edge did not present 
as much of a problem for Cycle 4. Figure 7 shows the kind of solution that can be 
expected for the same case as above. This figure was generated without any of the 
options discussed in the following paragraphs. 

In Cycle 4, the user may request the grid generator to use a double spline 
interpolation, one for each of the upper and lower surfaces, that will guarantee the 
midpoint in the C direction to be positioned at the leading edge. This feature 
was added because the flow solver makes an assumption about the flow direction at 
the wing surface when its magnitude is supersonic. It arose from the considerations 
that some wing cross sections have different upper and lower arc lengths, and that a 
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large velocity gradient is expected near a sharp leading edge, especially at nonzero 
angle of attack. To select this option in the grid generator, set SLEG =1.0 in 
the input list. 

Assuming that the above option has been taken, another problem arises. Call 
the leading edge point i = IHALF, the wing surface at k = NK, and assume that j 
ranges along wing span stations. With the consistent differencing scheme, three 
metrics at (IHALF, j+1/2, NK) must be computed. It is not unusual for such points 
to have grid lines leaving the surface at obtuse angles. Since one-sided differ- 
ences are still used at inner boundary points, the Jacobian may be numerically zero 
or even slightly negative at such points. Given a sophisticated grid generator, 
this problem can be eliminated. In the present case, however, a different approach 
was taken. The metrics and fluxes at the leading edge are split into upper and 
lower contributions, located at (IHALF-1/4, j+1/2, NK) and (IHALF+1/4, j+1/2, NK), 
respectively (identical to trailing edge treatment). A study involving the M6 wing, 
built into TWING, shows negligible impact on a wing with a blunt leading edge, while 
some improvement has been observed for a wing with a sharp leading edge. Finally, 
the user has the option of averaging the densities at the (IHALF- 1/2, j, NK) and 
(IHALF+1/2, j, NK) points with their two neighbors at IHALF-3/2 and IHALF+3/2. 

To turn this on, set SLET =1.0 in the input list. This leads to a more accurate 
plot of the coefficient of pressure in the vicinity of the leading edge for sharp 
leading edge cases. 


COMPARISONS WITH A TWO-DIMENSIONAL FLOW SOLVER 


The 2-D flow solver developed by Flores et al. in reference 1 was a modified 
version of TAIR, described in reference 5. It is a 2-D full-potential transonic 
flow code for airfoils, and is very similar to TWING. The features of the former 
flow solver are embedded in the code known at NASA Ames Research Center as 
TAIR10K. It was used in order to make two final comparisons. 

Using the same cross section as that of the wing presented in figures 6 and 7, 
we show the effect of higher Mach number. For these comparisons, the wing grid 
cross sections of TWING and the TAIR10K grid have the same number of points, 

151 x 31, and the wing has aspect ratio 20.0. Figure 8 shows a good comparison 
between the two codes for Mach 0.70. The flow on the upper surface is slightly 
supersonic without a shock. The wing code does not exhibit the peak at the leading 
edge as does the airfoil code. This is a consequence of the double interpolation at 
the leading edge lower and upper surfaces in TWING. Figure 9 shows what happens at 
Mach 0.75. In this case, a strong shock forms at about 80# of chord. The pressure 
peak at the leading edge for the 2-D flow solver is observed again. It is interest- 
ing that the shock for the wing is located slightly forward of the shock for the 
airfoil. It probably is caused by a 3-D relief effect in the wing calculation. 
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Figure 8.- Sharp leading edge, = 0.70, a = 0. 
( ) THING Cycle 4. ( ) TAIR10K. 



x/c x/c 

Figure 9.- Sharp leading edge, M,,, = 0.75, a = 0. 
( ) THING Cycle 4. ( ) TAIR10K. 

CONCLUSIONS 


The full-potential, steady, transonic, wing-flow solver, THING, has been modi- 
fied so that free-stream density and vanishing residual are captured in regions of 
constant velocity. Numerically precise consistency, which does not depend on the 
fineness of the computational grid, is obtained by slightly altering the differenc- 
ing scheme without affecting the implicit AF2 solution algorithm. 
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The computational efficiency remains about the same as before because the 
changes chiefly affect the fifteen metrics per grid point, which are computed once 
and stored. For each iteration, operations that were added for data expansion are 
offset by averaging operations that were eliminated. The number of iterations to 
convergence has been observed to increase for some cases, but decrease for others. 

With this new method, the outer boundary condition is accurately captured. We 
have seen that the smoothness of the solution is especially improved near regions of 
grid discontinuity. 
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